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We apply a Bayesian data analysis scheme known as the Markov Chain Monte Carlo (MCMC) 
to the tomographic reconstruction of quantum states. This method yields a vector, known as the 
Markov chain, which contains the full statistical information concerning all reconstruction param- 
eters including their statistical correlations with no a priori assumptions as to the form of the 
distribution from which it has been obtained. From this vector can be derived, e. g. the marginal 
distributions and uncertainties of all model parameters and also of other quantities such as the 
purity of the reconstructed state. We demonstrate the utility of this scheme by reconstructing 
the Wigner function of phase-diffused squeezed states. These states posses non-Gaussian statistics 
and therefore represent a non-trivial case of tomographic reconstruction. We compare our results 
to those obtained through pure maximum-likelihood and Fisher information approaches. 
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I. INTRODUCTION 



The tomographic reconstruction of quantum states 
represents an important laboratory tool in both quantum 
optics and quantum information alike. It can be applied, 
for example, to the reconstruction of the dynamical inter- 
action between quantum systems in a technique known 
as quantum process tomography. This latter technique is 
important in quantum computing, where the characteri- 
zation of quantum gates is essential to the overall quan- 
tum circuit Q, Hjj. Tomography has even been used 
as part of the optimization of advanced interferomctry 
as was done in the preparation of frequency dependent 
squeezing 

Since the theoretical discovery by Vogel and Risken 
[H that the Wigner function can be reconstructed from 
homodyne detector data, a number of reconstruction 
schemes have been developed ranging from direct inver- 
sion of the tomographic data by means of the filtered- 
back projection method @ to statistical methods such 
as maximum- likelihood estimation @, H, H, QjJ HH, [HI- 
An important feature of maximum-likelihood methods 
is the guaranteed positive "semi-definiteness" of the re- 
constructed state. The result of a maximum-likelihood 
reconstruction method is either a density matrix floL IT]| 
or a set of parameters (l3l [l4| which have maximized the 
likelihood functional given a model of the measurement 
apparatus and of the parameterized state. A full analysis 
of the experimental data, however, should also answer im- 
portant questions regarding error-bars on the estimation 
of the state parameters, possible correlations amongst the 
parameters, and error-propagation when using the recon- 
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structed state for further calculations of quantities such 
as the purity of the state or amount of entanglement. 

Several methods have been proposed in the literature 
to put error-bars on the reconstructed quantum states. 
For linear reconstruction techniques based on the aver- 
aging of sampling functions over the experimental data 
[15l | one can calculate the statistical uncertainties of the 
estimated quan tities by evaluating variances of the linear 
estimators [la ]. Another possible approach is to numer- 
ically simulate the whole measurement and reconstruc- 
tion process many times assuming that the reconstructed 
state is the true state that is being measured upon. The 
error bars are then calculated from the resulting ensem- 
ble of reconstructed states 17]. Finally, uncertainties 
on estimates obtained by maximum-likelihood method 
can be determined by evaluating the Fisher information 
matrix [l8[ . This latter approach essentially relies on ap- 
proximation of the likelihood function by a Gaussian and 
becomes exact only asymptotically in the limit of a very 
large number of experimental data. 

In the present paper we show that the uncertainties 
in quantum state estimation can be consistently deter- 
mined by using a general and statistically well motivated 
Bayesian analysis scheme known as the Markov Chain 
Monte Carlo (MCMC). The method is based on the im- 
plementation of a Markov chain to search the parameter 
space resulting in a set of samples from the joint poste- 
rior probability density distribution on the unknown pa- 
rameters of the model. This technique produces several 
important results. First, it yields the Markov chain con- 
taining all of the relevant statistical information about 
the parameter space. Second, one can extract a set of 
marginalized probability density distributions for each 
parameter quantifying the degree of uncertainty on their 
estimation. Third, the resulting chain can be used in 
further calculations, where one can produce probability 
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density distributions on quantities such as of the purity 
or amount of entanglement of the reconstructed state. 

This paper is divided into the following sections: in 
Sec. |n] the necessary concepts required from Bayesian 
data analysis are introduced and applied to the case of 
quantum state estimation. In Sec. IIIH the quantum like- 
lihood function for phase-diffused squeezed states is de- 
rived. In Sec. IIVI the Markov chain Monte Carlo algo- 
rithm is introduced and finally in Sec. [V] both the techni- 
cal details of the experimental realization as well as the 
results of the reconstruction are presented. 



II. BAYESIAN DATA ANALYSIS 

Baycs' Theorem prescribes the rule to invert the rela- 
tionship between the experimental data already observed 
and the parameterized model which could have generated 
the measured data set. The theorem reads 



and reads 
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where we use D to represent our data, / to represent 
our prior information or model of the experiment and A 
to represent a vector of model parameters. The quan- 
tity p{X\D ,1) is known as the posterior distribution, 
p(D\X,I) is known as the likelihood function, p(X\I) is 
the prior distribution and p(D\I) is a normalization fac- 
tor. The standard application of Bayesian analysis is to 
calculate the posterior distribution given a parameterized 
model of the sought after signal, the measured data and 
prior probability distributions on the values of the model 
parameters. These three elements are brought together 
through Eq. Q. 

Let us now construct the likelihood function for the 
quantum estimation problem. A general measurement 
on a quantum system can be described by the so-called 
positive operator- valued measure (POVM). Each possi- 
ble measurement outcome j is associated with a POVM 
element YLj which is a positive semidefinite operator. 
The probability of outcome j can be calculated as Pj = 
Tr[Hj p(X)] where p(A) denotes the density matrix of the 
measured quantum system that depends on the model 
parameters A. Since the total probability of some out- 
come is 1, the POVM elements sum up to identity oper- 
ator, Y^, j Hj = 1. This generic framework in particular 
encompasses a tomographic reconstruction of the state 
p(X) that consists of several different measurements M 
with possible outcomes indexed by Im- Then j = (M, Im) 
becomes a multi-index indicating both the measurement 
setting and the measurement outcome for a given setting. 
Let rij denote the observed number of measurement out- 
come j and N = . nj represents the total amount of over the bin, 
collected data. The likelihood function L is the probabil- 
ity of observation of a particular set {nj} for a given A. 
It follows that C is given by a multinomial distribution 



In terms of the constituents of Bayes's theorem the theo- 
retical probabilities {Pj} are functions of the parameters 
which are to be determined. The measured numbers of 
counts {nj} correspond to the data. 



III. THE LIKELIHOOD FUNCTION FOR 
PHASE-DIFFUSED SQUEEZED STATES 

The phase-diffused squeezed states were first intro- 
duced within the context of continuous-variable squeez- 
ing purification by [ijj HO, Hl[ . They arise when squeezed 
states are transmitted over de-phasing quantum channels 
such as optical fibers affected by thermal fluctuations. 
These states are characterized by non-Gaussian statistics 
and therefore represent a non-trivial case for quantum 
state tomography. 

The Wigner function for phase-diffused squeezed state 
is given by 
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with 



where x^ = x cos cf> + p sin </>, p^ = p cos 
x and p as the standard position and phase quadratures 
and 4> representing the random phase shifts distributed 
according to some probability distribution $ (0) . We use 
V x to represent the variance of the squeezed quadrature 
and V p for the variance of the anti-squeezed quadrature. 
We normalize the variances such that for vacuum state 
we have V x = V p = 1 and the state is squeezed in x 
quadrature if V x < 1. In the experiment we measure sev- 
eral different rotated quadratures xg, where 8 defines a 
specific measurement setting. The theoretical homodyne 
probability density distribution p(xg) can be calculated 
from Wigner function as a marginal distribution. In- 
tegration of W{x,p) over the conjugate quadrature pe 
yields, after some algebra 
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where V(4>) = V x cos 2 4> + V p sin 2 0. 

The data from each measurement is binned into L 
bins whose lower boundaries are defined by Qe.i- The 
outer bins extend to infinity and we set Qe,\ — — oo and 
Qe,L+i = °° ■ The corresponding theoretical probability 
Pe l is given by integration of the probability density Q 



Qe,i+i 



p(x e )dxg. 
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In our experiment, the results from two quadrature mea- 
surements were formed into histograms each containing a 
total of L = 70 bins. From the perspective of direct data 
inversion this corresponds to an overdetermined system, 
because we need to estimate only three real parameters, 
c. f. below. 

The POVM elements describing such binned homo- 
dyne detection can be expressed as 



\xg){xg\dxg, 



(6) 



where \xg) is an eigenstate of quadrature operator xg. 
Note that, by definition, the sum of theoretical probabil- 
ities over all bins is equal to one, 



(7) 



This is a mathematical expression of the fact that the ho- 
modyne detection always yields some outcome and, after 
each measurement, one of ng^i is increased by one. Put 
in a different way, the homodyne detection is described 
by a complete POVM ((6|) whose elements Hgj satisfy the 
condition J^i ^6,1 = 1- 

Assuming the phase noise distribution, $ (0), is a zero 
mean Gaussian, the state can be completely character- 
ized by just three parameters A = {V x , V p , V^} where V$ 
is the variance of the random phase shifts. The quantum 
log-likelihood function for the phase diffused squeezed 
states is finally obtained by taking the natural logarithm 
of Eq. @ giving us 



A = y^rtg,; In 



(8) 



where, for simplicity, we ignore all terms that do not 
depend on the parameter values. 



IV. MARKOV CHAIN MONTE CARLO 
ALGORITHM 
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FIG. 1: (Color online) Experimental setup: The experimen- 
tal setup consists of an optical parametric amplifier (OPA) for 
the generation of the squeezed vacuum states, a phase shifter 
to induce the random phase noise and a homodyne detector 
to measure the prepared state. The mode cleaner was used 
to increase the fringe contrast at the balanced homodyne de- 
tector. 



the choice of uniform priors has negligible effect on the 
numerical results of the MCMC [2{|. In this case Bayes's 
theorem, Eq. |T|), tells us that the likelihood function is 
proportional to the posterior distribution and since this 
is a function that we can compute for a given parameter 
space location we can use a standard sampling algorithm 
such as the Metropolis-Hastings sampler to draw from 
it. We describe our implementation of this sampler in 
Appendix fA] 



The goal of our Baycsian reconstruction scheme is the 
calculation of marginalized posterior distributions on our 
model parameters A = {V x , V p , V^}. To this end, the 
MCMC method can be used to generate samples drawn 
from the posterior distribution p(X\D,I). Since this dis- 
tribution is unknown one cannot directly sample from it 
and instead we sample from the distribution that is the 
product of the likelihood and the prior. The prior is cho- 
sen by considering the possible values of the parameters 
to be determined. Since the parameters to be determined 
in this case are variances, their values must be greater 
than zero. In order to assume relative ignorance in the 
value the parameters could take, we use a prior which 
only requires the variances to be positive and satisfy the 
Heisenbcrg uncertainty relation, V x V p > 1. Since the 
likelihood in this analysis is a sharply peaked function 



V. EXPERIMENTAL IMPLEMENTATION 
A. Description of the Experiment 

Figure [1] shows the experimental setup that was used 
to prepare the phase diffused squeezed states. The full 
details of the setup are provided in [2p| and will be sum- 
marized here. The squeezing source was an optical para- 
metric amplifier (OPA) constructed from a type I non- 
critically phase-matched MgO:LiNb0 3 crystal inside a 
standing wave resonator, similar to the design that previ- 
ously has been used in [22[ . The OPA was pumped with 
50 mW of green light at 532 nm resulting in a classical 
gain of about 1 1 . The length of the OPA cavity as well as 
the phase of the second harmonic pump beam were con- 
trolled using radio-frequency modulation/demodulation 
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FIG. 2: Markov chains: This figure depicts the evolution of the Markov chains where the abscissa represents the iteration number 
and the ordinate represents the value of the chain. After an initial "burn-in" period, in which the chains head to their steady- 
state positions, the chains eventually converge to a region of parameter space and begin to sample from the posterior distribution. 
The proposal distributions were taken to be Gaussian with the standard deviations = 0.0042, S p = 0.022, = 0.0037 
corresponding to the squeezing parameter, the anti-squeezing parameter and the phase noise parameter, respectively. The 
starting values were randomly chosen with the only constraint that they be non- negative and obey V X V P > 1. The chains 
settled to their equilibrium positions with means of fi x — 0.316, /j. p — 6.888, and /x^ = 0.171 corresponding to the squeezing, 
anti-squeezing and phase noise parameter respectively. 




FIG. 3: Marginalized posterior distributions: Each probability density was calculated using every tenth point from the marginal- 
ized chains corresponding to 38000 data points. The standard deviation of each density are a x = 0.0056, a v — 0.0289, and 
(70 = 0.0020 corresponding to the squeezing parameter, the anti-squeezing parameter and the phase noise parameter, respec- 
tively. The dotted line represents the result of the Fisher analysis. 



techniques. The mode cleaner was operated in high fi- 
nesse mode T = 10500 resulting in a line width of 55 kHz. 
A non-classical noise power reduction of slightly more 
than 5.0 dB was directly observed with a homodyne de- 
tector in combination with a spectrum analyzer at a 
Fourier sideband frequency of 6.4 MHz. 

The phase noise was induced by reflecting the squeezed 
field from a piezo-electric transducer (PZT) mounted 
high-reflection mirror that was quasi-randomly moved. 
The voltages applied to the PZTs were produced as fol- 
lows. An independent random number generator pro- 
duced data strings with a Gaussian distribution. The 
strings were digitally filtered to limit the frequency band 
to 2-2.5 kHz. The output interface was a common PC 
sound card with SNR of -HOdB. The sound volume was 
set to meet the desired standard deviation of channel 



phase noise. 

Homodyne detection confirmed that the squeezing 
degraded in the same way when phase noise was in- 
creased. The detector difference currents were electroni- 
cally mixed with a 6.4 MHz local oscillator. The demodu- 
lated signals were then filtered with a 400 kHz bandwidth 
low-pass filter and sampled with one million samples per 
second and 14 bit resolution using a National Instruments 
analog-digital sampling card. 



B. Results of the MCMC 

Figure [5] depicts the resulting chains after 40,000 it- 
erations of the MCMC algorithm. The abscissa repre- 
sents the number of iterations of the MCMC whereas the 
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ordinate represents the parameter values. After an ini- 
tial "burn-in" period of approximately 1000 iterations, 
in which the chain heads towards equilibrium, the chain 
converges and begins to sample from the posterior dis- 
tribution (which in this Gaussian case also includes the 
region of maximum-likelihood). The development of cri- 
teria for the determination of chain convergence is a gen- 
eral problem which has been the subject of much research 
[23l |24j • The general idea is to run multiple chains per 
estimation parameter and monitor their evolution both 
within each change and across each change. Convergence 
is inferred if all chains behave consistently. With respect 
to the case at hand, convergence of the chain can be in- 
ferred by comparing the locations to which the marginal- 
ized parameter chains have settled with the independent 
measurement of the V x , and V p parameter values per- 
formed with a spectrum analyzer. In the absence of such 
an independent measurement, the criteria in [23l. |24| can 
be used to infer convergence. 

The width of the marginalized chains, i. e., their stan- 
dard deviations, quantify the degree of uncertainty on 
the value of each parameter. By forming histograms of 
the chain as a function of each of the parameters we 
obtain their marginalized posterior probability distribu- 
tions as shown in Fig. [3] From these posteriors we ob- 
tain the following uncertainties on the model parameters: 
(T x = 0.0056 for the squeezing parameter, a p = 0.0289 
for the anti-squeezing parameter and finally er^ = 0.0020 
for the phase noise parameter. The proposal distribu- 
tions, from which the posterior distribution samples have 
been drawn, were chosen to be Gaussians with the fol- 
lowing standard deviations S x = 0.0042, S p = 0.022, 
= 0.0037 corresponding to the squeezing parameter, 
anti-squeezing parameter and phase noise parameter, re- 
spectively. These values were obtained through manual 
tuning of the MCMC algorithm. This is done by adjust- 
ing the individual standard deviations, i. e. S x , S p , S^, 
until the proportion of accepted jumps reaches approxi- 
mately 44% [25jj. 

As an independent test of the posterior standard devi- 
ations, we also calculated the Fisher information matrix 
IH, [27| given by 

_/dAdA\ 

The inverse Fisher matrix represents the covariancc of the 
posterior probability distribution for the true parameters 
A as inferred from a single experiment assuming Gaus- 
sian noise and constant priors over the parameter range 
of interest. The calculated standard deviations from the 
Fisher matrix are presented in Fig. [3] as well as in Table 
U where very good agreement is readily seen. It should 
be stressed that the posterior distribution standard de- 
viations obtained from the Fisher matrix are only valid 
when assuming Gaussian noise whereas the posterior dis- 
tribution standard deviations obtained from the Markov 
chain is valid regardless of the form of the posterior dis- 
tribution. The results of the MCMC as well as of the 



TABLE I: Standard deviations of posterior distributions: 
This table compares the standard deviations of the param- 
eter posterior distributions obtained from the Markov chain 
Monte Carlo method and from the Fisher information matrix. 
The standard deviations represent the error on the estimation 
of the parameter values. The Fisher matrix returns the actual 
standard deviations only in the case of Gaussian noise. 





Parameter Estimates 


Uncertainties 


Parameter 


MCMC Maximum Likelihood 


MCMC Fisher 


Vx 


0.316 0.317 


0.0056 0.0055 


v P 


6.889 6.880 


0.0289 0.0294 




0.171 0.171 


0.0020 0.0020 



Fisher analysis are compared in Table fl] 

Figure 2] depicts the evolution of the log-likelihood 
function (Eq. ©) for each iteration of the MCMC. It 
is seen that as the chains evolve through the "burn-in" 
stage the log-likelihood quickly increases. After approxi- 
mately 1000 iterations it has reached equilibrium where- 
upon parameter space jumps to higher likelihood values 
are balanced by jumps to lower likelihood values. 

Figure [5] represents the spectrum of the squeezed state 
before the phase diffusion. The measured state is seen 
to have a squeezing strength of -4.98 dB or a variance of 
V x = 0.31 and an anti-squeezing strength of 8.39 dB or a 
variance of V p = 6.91. These values have not been cor- 
rected for dark noise and were measured using a video 
bandwidth (VBW) of 10.0 Hz, a resolution bandwidth 
(RBW) of 100 kHz and a sweep time (SWT) of 1.5 s. 
These values lie within the width of the respective pos- 
terior distributions obtained from the MCMC analysis. 
Furthermore, only two quadrature measurements, each 
containing just 100, 000 data samples, were required to 
obtain these results. This represents a significant sav- 
ings in terms of experimental effort to reconstruct a non- 
Gaussian state. 

We also compare our results with that of a pure 
maximum-likelihood approach. This was achieved by al- 
tering step 3(d) of the Metroplis-Hastings sampler (see 
Appendix [XJ such that the chain remains at its cur- 
rent position whenever the likelihood ratio r is less than 
one. This forces the chain to only move to regions of 
higher likelihood and hence quickly converge on the true 
maximum-likelihood location. The parameter values ob- 
tained at the maximum-likelihood were: V x — 0.317 
for the squeezing parameter, V p — 6.880 for the anti- 
squeezing parameter and Vj, = 0.171 for the phase noise 
parameter. These are, however, the only results obtained 
from the maximum-likelihood estimation alone. Since 
the MCMC chain contains a complete statistical descrip- 
tion of the parameters, the statistical error on the re- 
construction of both the quantum state itself as well as 
derived quantities from it, e. g. purity, can be exactly 
determined. This will be illustrated in the next section. 
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C. Reconstruction of the Quantum State 



10 r 



Having completely characterized the parameter space, 
the final state can be reconstructed. In order to gen- 
erate the Wigner function shown in Fig. [HI the Markov 
chain together with Eq. ([3]) were used to calculate the av- 
erage Wigner function. The characteristic non-Gaussian 
statistics of the phase-space Wigner distribution is clearly 
manifested in the reconstructed state. It should be 
stressed that this averaging already takes into account 
the standard deviation of each value of the Wigner func- 
tion since the parameter values are taken from the pos- 
terior distribution. 

In addition to generating a phase-space plot of the 
reconstructed quantum state, the Markov chain can be 
used in further calculations of such properties as the pu- 
rity of the reconstructed state. Figure [7] depicts such a 
result. Using the analytical definition of the purity 



anti-squeezing 
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47r J J W 2 (x , p) dxdp, 



(10) 



and the resulting chain from the MCMC, the purity can 
be calculated, automatically taking into consideration 
the statistical error and correlation on the parameters 
determined by the MCMC. The result is a probability 
density whose standard deviation quantifies the degree of 
uncertainty on the estimation of the purity. For the state 
in question, we obtain a purity of n = 0.5649 ± 0.0028. 
It is important to note that this information is delivered 
directly from the MCMC itself; no additional assump- 
tions as to the distribution of the errors and their corre- 
lation properties need to be made. Additionally, any one- 
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FIG. 5: Zero span measurement of squeezing and anti- 
squeezing: The recorded amount of squeezing and anti- 
squeezing as measured by a balanced homodyne detector and 
spectrum analyzer at a Fourier frequency of 6.5 MHz. It is 
seen that approximately -4.98 dB of squeezing and 8.39 dB 
of anti-squeezing were directly measured without dark noise 
correction. These correspond to variances of V x — 0.31 and 
V p = 6.91, respectively. The units of the vacuum have been 
set to 1 corresponding to OdB. 
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FIG. 6: (Color Online) Reconstructed Wigner function: The 
Wigner function of the measured state as reconstructed using 
Eq. (|3} and the Markov chain. This reconstruction contains 
all the relevant statistical information. 



FIG. 4: Log-Likelihood: Evolution of the log-likelihood as a 
function of iteration. The abscissa represents the iteration 
number and the ordinate the log-likelihood value. After ap- 
proximately 1000 iterations the log-likelihood appears to have 
attained it's maximum value at which point the chains have 
reached equilibrium and are sampling from the posterior dis- 
tribution. 



dimensional quantity can be calculated in this manner. 
For example, if estimating the amount of entanglement 
of a non-Gaussian state, the logarithmic negativity [28| 
can be calculated over the span of the resulting chain. 
The result will be a probability density quantifying the 
uncertainty on its value. 
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FIG. 7: Reconstructed purity posterior distribution: The pu- 
rity was calculated using the resulting Markov chain of 38,000 
parameter values and Eq. (|10|l . This results in a probability 
density function which exactly quantifies the uncertainty on 
the estimation of the purity. The purity is calculated to be 
fj, = 0.5649 ± 0.0028. No such result is possible with a pure 
maximum-likelihood approach. 
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FIG. 8: (Color online) MCMC program flow chart: A 
schematic explanation of the process of selecting MCMC sam- 
ples from the posterior distribution. The lower loop iterates 
until the chain has converged to its equilibrium position and 
begins to sample from this posterior distribution. 



VI. CONCLUSION 

We have applied a Bayesian data analysis scheme 
known as Markov chain Monte Carlo (MCMC) to the 
tomographic reconstruction of quantum states. Taking 
phase diffused squeezed states as an example, we have 
provided the details as to the derivation of the likelihood 
function as well as to the numerical implementation of the 
MCMC. The results include a set of probability density 
distributions which exactly quantify the degree of uncer- 
tainty on the estimation of the parameters. These results 
were compared to both a pure maximum-likelihood and 
a Fisher information approach. Furthermore, using the 
Markov chain in the calculation of the state's purity en- 
abled the construction of a probability density distribu- 
tion on the value of the purity, thereby quantifying the 
degree of uncertainty on its calculation. 

We note that MCMC scheme is completely general and 
can be applied to higher dimensional problems, such as 
the reconstruction of the density matrix, and will be the 
topic of future publications. 
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APPENDIX A: METROPLIS-HASTINGS 
SAMPLER 

Our implementation of the MCMC is based on the 
Metroplis-Hastings sampler. It is defined by performing 
the following steps (see also Fig. [8]) : 

1. Generate initial values for parameters, A 

2. Compute the quantity £q = £(D\Xq) 

3. Iterate the following over the index t until the chain 
has converged 

(a) Generate a trial parameter vector £ according 
to a proposal distribution g(Af_i) 

(b) Compute the quantity £ t = C(D\£). If 
V X V P < 1 or some Vj < then set £ t = 0. 

(c) Compute the ratio r = £ t /£t-i 

(d) Sample from a uniform distribution, U(0 , 1) 
.r(r>U, set At = £, £t = £ accepted 

\ r < U, set A t = At_i, £t = £t-i rejected 
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The proposal distribution q(X) used in stage 3(a) is 
used to select trial parameter values within the MCMC. 
Theoretically q can be any distribution, however in prac- 
tice it is sensible for the proposal distribution to suggest 
jumps that are local to the current location but large 
enough to allow an efficient exploration of the parame- 
ter space [3(J. Stage 3(b) ensures that the sampling is 



restricted to subspace of physically admissible values of 
parameters A. The simple difference between this MCMC 
method compared to that of a pure maximum- likelihood 
method is that for repeated stages within the Metropolis- 
Hastings sampler the chain has finite probability of jump- 
ing both to higher or lower values of likelihood. 
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